/*
 * =====================================================================================
 *
 *       Filename:  getnonzeros.cc
 *
 *    Description:  Print out nonzero entries of either Q, A, or M operators (field index)
 *
 *        Version:  1.0
 *        Created:  06/17/2010 10:21:25 AM
 *       Revision:  none
 *       Compiler:  gcc
 *
 *         Author:  YOUR NAME (), 
 *        Company:  
 *
 * =====================================================================================
 */
#include <stdlib.h>
#include <vector>
#include <string>

using namespace std;

struct rank3 {
	char i; char j; char k;
	rank3(int e1, int e2, int e3){
		this->i = e1; this->j = e2; this->k = e3;
	}
	bool operator==(const rank3 & other) const{
		return other.i == this->i && other.j == this->j && other.k == this->k;
	}
};

struct rank2 {
	char i; char j; 
	rank2(int e1, int e2){
		this->i = e1; this->j = e2; 
	}
	bool operator==(const rank2 & other) const{
		return other.i == this->i && other.j == this->j ;
	}
};

string getQ(const rank2 & qindex,const string & kvar){
	static string twopisquared("TWOPI2");
	string nb = string("(-") + twopisquared + string("*") + kvar + string(")");
	if( qindex == rank2(1,1) ) return nb+" - F1";
	if( qindex == rank2(1,2) ) return "F1";
	if( qindex == rank2(2,1) ) return "F2";
	if( qindex == rank2(2,2) ) return nb+" - F2";
	if( qindex == rank2(3,3) ) return nb+" + F3";
	if( qindex == rank2(3,4) ) return "-F3";
	if( qindex == rank2(4,3) ) return "-F4";
	if( qindex == rank2(4,4) ) return nb+"-F4";
	if( qindex == rank2(5,5) ) return "1";
	if( qindex == rank2(6,6) ) return "1";
	return 0;
}

string getQI(const rank2 & qindex, const string & kvar){
	static string twopisquared("TWOPI2");
	string nb = string("(-") + twopisquared + string("*") + kvar + string(")");
	string b = string("(") + twopisquared + string("*") + kvar + string(")");
	string nA = "( -1.0/ ( " + b + string("*(") + b + string("+F1+F2) ) )");
	string  B = "(  1.0/ ( " + b + string("*(") + b + string("+F3-F4) ) )");
	if( qindex == rank2(1,1) ) return nA+"*(" + b + string(" + F2)");
	if( qindex == rank2(1,2) ) return nA+"*F1";
	if( qindex == rank2(2,1) ) return nA+"*F2";
	if( qindex == rank2(2,2) ) return nA+"*("+b+"+ F1)";
	if( qindex == rank2(3,3) ) return B+"*("+nb+" + F4)";
	if( qindex == rank2(3,4) ) return B+"*(F3)";
	if( qindex == rank2(4,3) ) return B+"*F4";
	if( qindex == rank2(4,4) ) return B+"*("+nb+"+F3)";
	if( qindex == rank2(5,5) ) return "1";
	if( qindex == rank2(6,6) ) return "1";
}


const vector<rank2> & getnonzeroQs(){
	static vector<rank2> list;
	if(list.size() == 0){
		list.push_back(rank2(1,1));
		list.push_back(rank2(1,2));
		list.push_back(rank2(2,1));
		list.push_back(rank2(2,2));
		list.push_back(rank2(3,3));
		list.push_back(rank2(3,4));
		list.push_back(rank2(4,3));
		list.push_back(rank2(4,4));
		list.push_back(rank2(5,5));
		list.push_back(rank2(6,6));
	}
	return list;
};

const vector<rank3> & getnonzeroPs(){
	static vector<rank3> list;
	if(list.size() == 0){
		list.push_back(rank3(1,1,1)); 
		list.push_back(rank3(2,2,2)); 
		list.push_back(rank3(3,3,3)); 
		list.push_back(rank3(4,4,4)); 
		list.push_back(rank3(4,6,4)); 
		list.push_back(rank3(5,5,3)); 
	}
	return list;
}

string convertint2string(int a){
	static char * temp = NULL;
	if(!temp){
		temp = new char[5000];
	}
	sprintf(temp, "%i", a);
	string rvalue(temp);
	return rvalue;
};

template <class rank>
bool doesindexexist(const rank & index, const vector<rank> & list){
	for(size_t i = 0; i < list.size(); i++){
		if( list[i] == index ) return true;
	}
	return false;
};


int printcppplus(){
	const vector<rank2> Qs = getnonzeroQs();
	const vector<rank3> Ps = getnonzeroPs();
	int num = 0;
	printf("switch(aiter){");
	for(int i = 1; i <= 6; i++){
		for(int j = 1; j <= 6; j++){
			for(int k = 1; k <= 6; k++){
				// internal sums
				string expr;
				for(int m = 1; m <= 6; m++){
					for(int n = 1; n <= 6; n++){
						bool nonzero = false;
						rank2 QIindex(i,m);
						rank2 Qindex1(n,j);
						rank2 Qindex2(n,k);
						rank3 Pindex1(m,n,k);
						rank3 Pindex2(m,n,j);
						if( doesindexexist(QIindex, Qs) ){
							string texpr = string("(") + getQI( rank2(i,m) , string("k2" ) ) + string(")");
							//string("Qinv[ rank2(iQ") + convertint2string(i) + convertint2string(m) + string(", ki) ]");
							string lexpr;
							string rexpr;
							if(doesindexexist(Pindex1, Ps) && doesindexexist(Qindex1, Qs)){
								nonzero = true; 
								lexpr = string("(") + getQ( rank2(n, j), string("p2") ) + string(")");
								//string("Q[rank2(iQ") + convertint2string(n) + convertint2string(j) + string(", pi)]");
							}
							if(doesindexexist(Pindex2, Ps) && doesindexexist(Qindex2, Qs)){
								nonzero = true; 
								rexpr = string("(") + getQ( rank2(n,k), string("q2") ) + string(")");
								//string("Q[rank2(iQ") + convertint2string(n) + convertint2string(k) + string(",qi)]");
							}
							if(nonzero){
								texpr += string("*(");
								if(lexpr.size() > 0) texpr += lexpr;
								if(rexpr.size() > 0) texpr += string("-") + rexpr;
								texpr += string(")"); 
								if(expr.size() > 0){
									expr += string("\n\t\t+") + texpr;
								}else{
									expr += texpr; 
								}
							}
						}
					}
				}
				if(expr.size() > 0){
					//expr = string("\\frac{1}{2}\\left\\{") + expr + string("\\right\\}");
					//if(num > 0) printf("\\\\");
					printf("\n\tcase iA%i%i%i:\n\t\treturn %s;", i,j,k, expr.c_str());
					num++;
				}
			}}}
	//printf("\\\\\n\\textrm{count} & = & %i\n", num);
	printf("\n};");
	return 0;
};

int printlatex(){
	const vector<rank2> Qs = getnonzeroQs();
	const vector<rank3> Ps = getnonzeroPs();
	int num = 0;
	printf("\\documentclass{article}\n");
	printf("\\begin{document}\n");
	printf("\\begin{eqnarray}");
	for(int i = 1; i <= 6; i++){
		for(int j = 1; j <= 6; j++){
			for(int k = 1; k <= 6; k++){
				// internal sums
				string expr;
				for(int m = 1; m <= 6; m++){
					for(int n = 1; n <= 6; n++){
						bool nonzero = false;
						rank2 QIindex(i,m);
						rank2 Qindex1(n,j);
						rank2 Qindex2(n,k);
						rank3 Pindex1(m,n,k);
						rank3 Pindex2(m,n,j);
						if( doesindexexist(QIindex, Qs) ){
							string texpr = string("Q^{-1}_{") + convertint2string(i) + convertint2string(m) + string("}(k^{2})");
							string lexpr;
							string rexpr;
							if(doesindexexist(Pindex1, Ps) && doesindexexist(Qindex1, Qs)){
								nonzero = true; 
								lexpr = string("Q_{") + convertint2string(n) + convertint2string(j) + string("}(p^{2})");
							}
							if(doesindexexist(Pindex2, Ps) && doesindexexist(Qindex2, Qs)){
								nonzero = true; 
								rexpr = string("Q_{") + convertint2string(n) + convertint2string(k) + string("}(q^{2})");
							}
							if(nonzero){
								texpr += string("\\left[");
								if(lexpr.size() > 0) texpr += lexpr;
								if(rexpr.size() > 0) texpr += string("-") + rexpr;
								texpr += string("\\right]"); 
								if(expr.size() > 0){
									expr += string("\n\t\t+") + texpr;
								}else{
									expr += texpr; 
								}
							}
						}
					}
				}
				if(expr.size() > 0){
					expr = string("\\frac{1}{2}\\left\\{") + expr + string("\\right\\}");
					if(num > 0) printf("\\\\");
					printf("\nA_{%i%i%i} & = & %s", i,j,k, expr.c_str());
					num++;
				}
			}}}
	//printf("\\\\\n\\textrm{count} & = & %i\n", num);
	printf("\\end{eqnarray}\n\\end{document}\n");
	return 0;
};

int printcpp(){
	const vector<rank2> Qs = getnonzeroQs();
	const vector<rank3> Ps = getnonzeroPs();
	int num = 0;
	printf("switch(aiter){");
	for(int i = 1; i <= 6; i++){
		for(int j = 1; j <= 6; j++){
			for(int k = 1; k <= 6; k++){
				// internal sums
				string expr;
				for(int m = 1; m <= 6; m++){
					for(int n = 1; n <= 6; n++){
						bool nonzero = false;
						rank2 QIindex(i,m);
						rank2 Qindex1(n,j);
						rank2 Qindex2(n,k);
						rank3 Pindex1(m,n,k);
						rank3 Pindex2(m,n,j);
						if( doesindexexist(QIindex, Qs) ){
							string texpr = string("Qinv[ rank2(iQ") + convertint2string(i) + convertint2string(m) + string(", ki) ]");
							string lexpr;
							string rexpr;
							if(doesindexexist(Pindex1, Ps) && doesindexexist(Qindex1, Qs)){
								nonzero = true; 
								lexpr = string("Q[rank2(iQ") + convertint2string(n) + convertint2string(j) + string(", pi)]");
							}
							if(doesindexexist(Pindex2, Ps) && doesindexexist(Qindex2, Qs)){
								nonzero = true; 
								rexpr = string("Q[rank2(iQ") + convertint2string(n) + convertint2string(k) + string(",qi)]");
							}
							if(nonzero){
								texpr += string("*(");
								if(lexpr.size() > 0) texpr += lexpr;
								if(rexpr.size() > 0) texpr += string("-") + rexpr;
								texpr += string(")"); 
								if(expr.size() > 0){
									expr += string("\n\t\t+") + texpr;
								}else{
									expr += texpr; 
								}
							}
						}
					}
				}
				if(expr.size() > 0){
					//expr = string("\\frac{1}{2}\\left\\{") + expr + string("\\right\\}");
					//if(num > 0) printf("\\\\");
					printf("\n\tcase iA%i%i%i:\n\t\treturn %s;", i,j,k, expr.c_str());
					num++;
				}
			}}}
	//printf("\\\\\n\\textrm{count} & = & %i\n", num);
	printf("\n};");
	return 0;
};


int main(const int argc, const char * argv[]){
	for(int i = 0; i < argc; i++){
		if(string(argv[i]) == string("latex")){
			printlatex();
		}else if(string(argv[i]) == string("cpp+")){
			printcppplus();
		}else if(string(argv[i]) == string("cpp")){
			printcpp();
		}
	}
	return 0;
};
